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ABSTRACT 

We describe the dynamics of a stream of equally spaced macroscopic particles in orbit 
around a central body (e.g. a planet or star). A co-orbital configuration of small bod- 
ies may be subject to gravitational instability, which takes the system to a spreading 
disordered and collisional state. We detail the linear instability's mathematical and 
physical features using the shearing sheet model and subsequently track its nonlinear 
evolution with local N-body simulations. This model provides a convenient tool with 
which to understand the gravitational and collisional dynamics of narrow belts, such 
as Saturn's F-ring and the streams of material wrenched from tidally disrupted bodies. 
In particular, we study the tendency of these systems to form long-lived particle ag- 
gregates. Finally, we uncover an unexpected connection between the linear dynamics 
of the gravitational instability and the magnetorotational instability. 
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This paper is concerned with the gravitational and col- 
lisional dynamics of belts of material orbiting planets or 
stars. Of particular interest are Saturn's F-ring and the tidal 
streams torn from disrupted satellites. The F-ring is thought 
to comprise a population of large objects (of some 10 km) 
swathed in dust (Showalter 2004, Esposito et al. 2008, Mur- 
ray et al. 2008). Being located so near the Roche limit, the 
size distribution of these larger bodies evolves according to 
gravitational aggregation, and tidal and collisional disrup- 
tion (Barbara and Esposito 2002, Esposito et al. 2012). The 
tidal environment is very different for those dense narrow 
rings located interior the Roche limit, such as the e-ring of 
Uranus or the dense ringlets ensconced in Saturn's C-ring 
(Colwell et al. 2009). It is possible that these dense ringlets 
originated from the tidal disruption of small satellites (or 
the mantles thereof), with their early evolution controlled 
by gravitational instabilities (Leinhardt et al. 2012). 

Wc would like to understand these several processes — 
tidal disruption, particle clumping, gravitational instability, 
collisions — in a single theoretical framework using a simple 
but illuminating analytical model and N-body simulations. 
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The main configuration that we consider is the gravitational 
stability of a stream, or line, of many equally spaced par- 
ticles. To facilitate its study, we employ a local model: the 
shearing sheet. In doing so we can more easily reproduce 
the famous result of Maxwell's Adams Prize essay (Maxwell 
1859, Cook and Franklin 1964) and form a more intuitive un- 
derstanding of the instability's physics. We find a stream of 
particles is susceptible to two kinds of disruption: a familiar 
gravitational clumping instability, and a growing epicyclic 
mode. Which instability the system selects depends on a 
single dimensionless parameter analogous to the Roche pa- 
rameter that controls the tidal disruption of a satellite. 

The nonlinear outcome of these instabilities is studied 
with N-body simulations using the code REBOUND (Rein 
and Liu 2012). We find that the precise form of the ini- 
tial instability is unimportant; instead, the nonlinear evo- 
lution is entirely controlled by the tendency of the system 
to form particle aggregates. Existing clustering criteria are 
discussed (Weidenschilling et al. 1984, Ohtsuki 1993, Canup 
and Esposito 1995), and we simulate an interesting interme- 
diate regime relevant to Saturn's outer rings, whereby col- 
lisional aggregation is rare, but not impossible, and instead 
particles form temporary short-lived clusters (cf. 'dynamical 
ephemeral bodies', Weidenschilling et al. 1984). 

These gravitational instabilities may illuminate other 



astrophysical contexts. For instance, the mechanism of in- 
stabihty we outline here could also be at work in the non- 
axisymmetric instabilities of narrow fluid rings and slen- 
der tori (Maxwell 1859, Goodman and Narayan 1988, Pa- 
paloizou and Lin 1989). It also manifests in the 'propeller 
and frog' resonance of Pan and Chiang (2010, 2012) which 
aims to describe the observed migration of large embedded 
objects (propellers) in Saturn's A-ring. If the propeller is 
permitted to back-react on the rest of the ring, then an insta- 
bility will arise, analogous to the one we study here. Lastly 
we show how the mechanism of instability in the growing 
epicyclic mode, relying on significant angular momentum 
exchange, shares some of the mathematical and physical 
characteristics of the magnetorotational instability (MRI) 
(Balbus and Hawley 1991). In fact, the two instabilities are 
nearly identical for a certain (somewhat artificial) equilib- 
rium set-up. 

The following section outlines the linear theory of the 
gravitational instability, connecting it to previous work by 
Maxwell (1859) and Fermi and Chandrasekhar (1953), while 
showing its similarities to the MRlM Section 3 presents N- 
body simulations of the instability's nonlinear and collisional 
evolution. In Section 4 we summarise our results. 



2 LINEAR THEORY 

2.1 Governing equations and setup 

Consider a train of spherical particles in a circular orbit 
around a point mass or around a central mass with a point- 
mass potential. The particles initially inhabit the same or- 
bital radius but are equally spaced in azimuth by a distance 
h. Each particle has the same mass m which is much smaller 
than the central mass. If the spacing h is sufficiently small 
and the ensuing dynamics remain confined to lengthscales 
much less than the orbital radius, we may adopt the shearing 
sheet model or local approximation (Goldreich and Lynden- 
Bell 1965), in which case we can write down the equations 
of motion for each particle in a convenient way. The shear- 
ing sheet is anchored at a radius _Ro from the central object 
around which it orbits with frequency Q. In this model, the 
motion of particle n is governed by 



Xn — 2Qy„ = 30, x„ + f", 
2/„ -|-2rii„ = fy, 



(1) 
(2) 
(3) 



where {x„, yn, z„) denotes the coordinates of particle n in 
the shearing box relative to the reference position Ro, an 
overdot indicates a time derivative, and f" is the specific 
collective gravitational force of all the other particles on par- 
ticle n. The gravitational force may be written as 



r=Gm^^^^"''^ 



a^n 



[Aa;|„-fAy2 -fAz?]3/2 



(4) 



^ At the conclusion of this work we discovered that some of the 
linear theory we present in Section 2 was treated in a similar way 
by Willerding (1986), though our emphases and interpretations 
differ. 



where the sum is over all other particles, with G as the 
gravitational constant, and the distances 



^Xjn 



^Vjn = yj -yn, 



In reality there are a large but finite number of particles in 
the ring, but in practice we take the above sum to be infi- 
nite and so j runs between — oo and oo. Even if the force is 
approximated as an infinite sum, it remains convergent be- 
cause only nearby particles significantly infiuence any given 
particle's dynamics. See Salo and Yoder (1988) or Vanderbei 
and Koleman (2007) for a discussion of situations in which 
the rings contain only a few massive particles. 



2.2 Equilibrium 

This set of equations admits as an equilibrium state a pro- 
cession of equally spaced particles located in the central axis 
of the sheet, i.e. for 



0, 



yn 



hn. 



0, 



VnG 



(5) 



where h is the azimuthal spacing. The force f" is zero for 
each n: the gravitational infiuences of the particles preceding 
a given particle n are cancelled exactly by all the particles 
trailing particle n. 



2.3 Linearised equations 

We next perturb this stream of particles by a set of small 
displacements {x'„, y'„) that lie only in the orbital plane. 
Vertical displacements are neglected as they decouple in the 
linear regime and do not lead to instability. The linearised 
equations for these time-dependent displacements are: 



Xji A\Ly^i o\L Xji 

y'n + 2fiiti = - 



Gm 



E 



^Gm y^ j/j 



yn 



ft3 



\j - n|3 



(6) 
(7) 



where both sums omit the j = n term. 

We now suppose the particle motions assume a collec- 
tive Fourier mode of the type 

/ Y" st'{-nhki I -\^ st-^-nhki 

x^ ^ e , y^ J e , 

with X and Y complex amplitudes, s a growth rate (complex 
in general), and k a (real) wavenumber. In order to avoid 
aliasing, the magnitude of k is restricted to be equal or less 
than n/h. The equations reduce to two, 

s'^x -2nsY -sn'^x = -^FX, (8) 



s'^Y + 2QsX = 2^FY, 



(9) 



in which we have introduced the dimensionless force function 
F. It can be manipulated into 



F{K) = 2J2 
3=1 



1 — cos(j K) 



(10) 



where we have introduced a dimensionless wavenumber K 



hk. The force function depends only on K, and can be re- 
expressed in terms of the zeta and Clausen functions (Lewin 
1981). We find that _F > and is an increasing function 
of K, with F — >■ as AT — )■ 0. The following asymptotic 
expansion, for small K, offers a reasonable approximation 
for all physical wavenumbers: 



F - 



^-InK 



K' 



96 



-0(K'^ 



(for a derivation see, for example, Nicorovici et al. 1994). 



2.4 Springs, clumping, and angular momentum 
exchange 

Before moving on to the dispersion relation itself, we make 
some preliminary remarks. First, if the central object and 
the orbital motion are neglected, the terms involving Q. in 
Eqs ([8| and S disappear. Consequently, the azimuthal and 
radial motions decouple and the azimuthal equation gives 
rise to an unstable longitudinal mode which tends to clump 
particles together. Its growth rate is clear from inspection 
and is equal to ^IGniF jK^ . This unstable mode, which 
grows as the result of a release of gravitational potential 
energy, is the discrete analogue of the instability of an in- 
finite self-gravitating cylinder described by Chandraskehar 
and Fermi (1953) (see also Chandrasekhar 1981). 

Second, (|8| and S are reminiscent of the equations 
governing two orbiting bodies attached by a spring, as pop- 
ularised by Balbus in his interpretation of the magnetorota- 
tional instability (Balbus and Hawley 1991). There is, how- 
ever, a crucial difference in that our 'gravitational spring 
constant' differs in both size and sign depending on the di- 
rection of the spring force. While the gravitational 'spring' 
seeks to resist extension in the radial direction (as in the 
MRI), it actively encourages it in the longitudinal direc- 
tion. In some sense, the spring wants to become 'clumpy'. 
At the same time the spring forces also facilitate angular 
momentum exchange, which introduces a separate route to 
instability, employed by the MRI and the Papaloizou-Pringle 
instability (Papaloizou and Pringle 1984). Outward angular 
momentum transfer liberates orbital energy which is chan- 
neled into exponential growth. The gravitational instabili- 
ties we study in this paper are the outcome of the interplay, 
and sometimes competition, between these two instability 
mechanisms. 



2.5 Dispersion relation 

We rescale time so that il — 1. Eliminating X and Y from 
Eqs (|8| and H yields the following dispersion relation for 
s: 

s* + [l- gF{K)]s^ + 2gF{K)[3 - gF{K)] = 0, (11) 

where we have introduced the important dimensionless 
quantity 



{m/h^) G 



(12) 



It measures the influence of the stream's self-gravity relative 
to its inertial forces. 



The g parameter bears a superficial resemblance to the 
'Roche parameter' of an orbiting fluid satellite vulnerable to 
tidal disruption: 



7^: 



_pG 



(13) 



where p is the mass density of the satellite. We find grav- 
itational instability favours larger g — hence greater mass 
densities and greater radii, Rq. In contrast, a larger 7?. in 
the Roche problem corresponds to a greater resistance to 
disruption (Chandrasekhar 1987). This reflects the fact that 
the inertial forces are stabilising in the ring context and dis- 
ruptive in the satellite context. As we show later, however, 
the parameter g does not control the nonlinear saturation 
of gravitational instability. 

Finally, note that by making the identification gF{K) — 



\k /Q, , where va denotes Alfven speed, Eq. (Ill be- 



comes strikingly similar to the MRI dispersion relation (e.g. 
Eq. (79) in Balbus 2003). But differences in sign occur at key 
places, which we associate with the tendency to azimuthally 
clump. 



2.6 Fastest growing modes and stability criterion 

In this problem the shortest modes are the most unstable: 
they always grow the fastest and they are always the flrst 
to be destabilised. Therefore, our stability considerations 
need only make reference to these short scales. The shortest 
modes possess K — n and the force function is, consequently, 
F{tv) = (7/2)C(3) « 4.2, where C(3) is Apery's constant. 
From the explicit solution to the (bi-) quadratic (111, we 
find three bifurcations as g varies. In order of increasing g, 
these occur at 



5 = 



26±8V10 
63C(3) 



and 



7C(3) 



0.00927, 0.677, 



0.714. 



(14) 



(15) 



The first bifurcation corresponds to the onset of insta- 
bility. When g < 0.00927 there exist only neutral modes, 
and the system is stable. Two of the modes in this regime 
consist of near-epicyclic motion with frequency |sj < SI. As 
g is increased through 0.00927, there is a Hopf bifurcation 
and both these oscillations pick up equal and positive growth 
rates (in addition, there are decaying modes as the system is 
Hamiltonian) . The gravitational instability takes the rather 
novel form of a pair of growing epicycles (as in the viscous 
overstability — Borderies et al. 1985, Schmit and Tschar- 
nuter 1995, Latter and Ogilvie 2009). Note that the stability 
criterion, g < 0.00927, can be reformulated into the follow- 
ing: 



^ < 2.298657/iV^ 



(16) 



where A'^ is the total number of particles and M is the mass 
of the central object. Here we have defined the number of 
particles through Nh = 2ttRo. Equation (16 1 is the famous 



stability criterion that Maxwell derived in his Adams Prize 



essay for the stability of A^ co-orbital particles for large A^ 
(Maxwell 1859, Cook and Franklin 1964). 

Now, when g is increased towards 0.677, the next critical 
value, the oscillation frequencies of the two unstable modes 
are gradually suppressed by the self-gravity. At g = 0.677 
these oscillation frequencies are precisely zero, and there ex- 
ists a double root to (111. When g > 0.677 both unstable 



modes are monotonically growing, and cease to oscillate. As 
g is increased further, these two modes decouple: one grow- 
ing more quickly than the other. Eventually, when g > 0.714 
(the third bifurcation), the slower of the two becomes neu- 
tral, leaving only one growing mode. The instability now 
takes the more familiar form of longitudinal clumping. This 
is especially clear in the limit of large g, where the unsta- 
ble mode's growth rate approaches ^^/2g F{'k), which is the 
value for a non-orbiting line of self-gravitating particles (cf. 
Section 2.4) . 




2.7 Modes of general wavelength 

Modes on longer scales (smaller K) follow the same pattern, 
though the critical values of g are larger, and are different 
for each K. Note that sufficiently long modes are always 
stable for any given g. Once we stipulate g — no matter 
how large — formally we can always find a sufficiently small 
K so that gF{K) ^ 1. (In reality, of course, K is limited by 



the circumference of the ring.) Expanding (111 in this small 
combination gives the growth rates: 



s = ±i 



[y/6gF + 0{gF)] , s ^ ±i[l + 0{gF)] . (17) 



The instability scale is tied to the strength of self-gravity. 
Conversely, in the MRI context, once the strength of the 
imposed magnetic field is specified (via the Alfven speed), 
we may always find a wavelength above which the MRI op- 
erates. In the gravitational instability, once the strength of 
the particles' self-gravity is specified, we may always find 
a wavelength below which the GI operates (but only when 
g > 0.00927). 



Figure 1. Growth rates of the gravitational instability as a func- 
tion of dimensionless wavenumber K for three different g = 0.01, 
0.08, and 0.25. Each curve is labelled with its corresponding value 
of g. There are two unstable modes with the same growth rate 
and both oscillating with a frequency S! itO. 




2.8 Growth rates 

In Fig. 1, we plot the real part of the growth rate s of the 
unstable modes as a function of K, for various g. In this 
figure we concentrate on the growing epicycles; thus we limit 
g to be less than g < 0.677. Three values of g are chosen: 
0.01, 0.08, and 0.25. As explained earlier, the shortest modes 
(largest K) modes are the most unstable, but K ^ n. Note 
finally that for sufficiently long modes instability is quenched 
in all three cases. 

In Fig. 2, we extend the previous plot to larger g, in 
which we find the instability taking the form of gravitational 
clumping on shorter scales. Two values oi g are chosen above 
the value 0.677. In both cases instability is quenched for 
small K. One may also observe the bifurcation at a second 
critical K at which point the instability changes from being 
two growing epicycles, to two monotonic clumping modes. 
For K greater than a slightly larger value, one of these modes 
stabilises, and only one growing mode remains. 



Figure 2. Growth rates as functions of K for larger values of 
g. Two values are chosen: g = 1, 10. Above a critical K, in 
both cases, there exists only one monotonically growing mode. 
For K below a second critical value, there are two unstable 
epicyclic modes growing at the same rate. This value occurs at 
the 'branches' in the dispersion relation. For K less than a third 
critical value the modes are stable. 



2.9 Physical interpretation 

The instability associated with longitudinal clumping, at 
larger g, is relatively straightforward to understand. A pat- 
tern of compression and rarefaction is self-amplifying be- 
cause gravitational attraction increases when neighbouring 
particles are displaced closer to each other. In contrast, the 
'epicyclic' instability, which occurs for smaller values of g, 
relies on a mechanism that is a little more subtle and which 
we now describe using a simple cartoon. 

Consider the four 'snapshots' of Fig. 3. In the first panel. 



i 



? 



particle, when displaced, feels the gravitational attraction 
of the other two, but does not act back upon them. Con- 
sequently, it undergoes a periodic epicyclic-like motion of a 
steady amplitude. However, as soon as we let the third par- 
ticle react back on its two chaperones, an instability of the 
type we describe above should occur. Recently this seems to 
have been observed (Pan and Chiang 2012). 



o 






o 




Figure 3. A cartoon of the onset of oscillatory instability. The 
first panel shows two adjacent particles given a small displace- 
ment toward one another. In the second, their mutual attraction 
gives rise to an epicyclic motion (akin to a 'horseshoe' turn), the 
amplitude of which is increased by a gravitational encounter with 
the particles on either side of the pair, shown in the third panel. 
In the last panel, the two particles kick each other once again 
giving rise to an even large (modified) epicyclic motion. 



we have drawn just two particles within the initial equilib- 
rium state. These two particles have been displaced slightly 
towards each other (and away from their other immediate 
neighbours). As a result, there is a net gravitational at- 
traction between them (the black arrows). This force 'slows 
down' the leading particle and 'speeds up' the trailing parti- 
cle; thus an angular momentum exchange takes place. Con- 
sequently, both particles begin an epicycle around the two 
new radii that are associated with their new angular mo- 
menta: the leading particle moves to a smaller radius, while 
the trailing moves to a larger radius (Panel 2). Note that 
if the gravitational attraction is sufficiently strong (large g) 
then the particles would clump or collide before the epicyclic 
motion gets properly underway. Panel 3 describes what hap- 
pens when the two particles complete half of their epicycle. 
At this stage they are closer to their other near neighbours 
(in white) than to each other. Both subsequently receive a 
second attractive force from the white particles that sends 
them into new (larger) epicycles. Finally, when our two black 
particles return to their initial radius (panel 4), they are az- 
imuthally closer to each other than when they began (in 
panel 1). Their gravitational attraction is now stronger and 
sends them into new epicycles of greater amplitude, and the 
process runs away. Thus the instability relies on a sequence 
of gravitational encounters between neighbouring elements; 
and these act as an 'epicyclic amplifier'. 

The amplifier we describe here should be contrasted 
with the 'propeller and frog' resonance (Pan and Chiang 
2010). In the set-up of Pan and Chiang, two 'particles' are 
held fixed ahead of and behind a third particle. This third 



2.10 A gravitational analogue of the MRI 

Before we move on to the nonlinear development of the grav- 
itational instability, it is of interest to briefiy examine a dif- 
ferent but related model of instability, in which the dynamics 
of the incompressible MRI are reproduced in most of their 
details. 

Consider, instead of a stream of particles extended 
along the azimuthal direction, a vertical line of particles lo- 
cated at a single radius and azimuth: Xn — i/n ~ 0. The 
particles are equally spaced in z. However, in order to keep 
this configuration in equilibrium it is necessary to 'switch- 
off' the vertical tidal force of the central object in Eq. ([3|. It 
is also necessary to assume that the vertical line of particles 
is extremely long and that we study only particles at the 
midplane far from the unbalanced ends. This 'equilibrium' 
may then be written as 



= 0, y„ = 0, 



= ftn. 



Vn e z. 



As before, we consider only linear planar disturbances to 
this basic state, and we describe this perturbation by 



(x'r,, y'„,0) = {X, Y, 0) 



st-\-nKi 



Its linearised equations read 



Gm 



FX, 



s^Y + 2QsX = ~^FY, 



(18) 
(19) 



where F is given earlier in Eq. ( 10 1 



Equations (18l-(19l describe two bodies connected by 
a spring, with a single spring constant {Gm/h?)F{K). The 
linear dynamics of the incompressible MRI can be described 
in precisely the same way if v\k^ is substituted for the spring 
constant (Balbus 2003). Under this substitution, the disper- 
sion relations of the MRI and gravitational instability are 
almost identical! The only difference comes from F's depen- 
dence on wavenumber, which is more complicated than k^ . 
Yet like k^ , it is an increasing function and goes to as 
fc^O. 

So in this toy problem we have managed to reproduce 
the MRI dynamics by letting the gravitational force between 
particles do the job of magnetic tension in the MRI. And 
by extending the particles in the z-direction but suppress- 
ing their vertical motion we can nullify gravity's tendency 
to clump the particles (the complicating process in Section 
2.4). Particles at different z are drawn apart from each other 
in the plane; but these displacements engender planar grav- 
itational torques that transport angular momentum from 
particles at smaller radii to particles at larger radii. Parti- 
cles then drift further apart, angular momentum is further 
exchanged, and the process runs away. The system develops 



a vertical sequence of planar jets, or 'channel flows', as in 
the MRI (Goodman and Xu 1984). Unhke the MRI, how- 
ever, these flows are not nonlinear solutions to the equations 
of motion. 



3 NONLINEAR AND COLLISIONAL 
EVOLUTION 

In this section we present three-dimensional N-body simula- 
tions of the nonlinear and coUisional evolution of the grav- 
itational instabilities discussed. The equations (|l])-(|3| are 
numerically evolved forward in time in a corotating local 
frame, which imposes a finite periodicity in the azimuthal 
[y) direction, but has no radial boundaries. Unlike the shear- 
ing box, particles cannot leave the box through one radial 
boundary and reappear through the other. The simulations 
can show how the saturation of the instability depends on 
the governing parameters of the system, and if the eventual 
outcome of the clumping form of the instability differs from 
the outcome of its oscillatory form. 



3.1 Parameters 

In addition to the dimensionless parameter g, which governs 
the linear stage, N-body simulations introduce three more 
parameters associated with the coUisional energy losses, the 
particle diameter, and the number of particles A'^ in the box. 
We, however, have verified that our results do not depend 
on A'^, which has no direct physical meaning. Particles are 
assumed to be non-spinning inelastic spheres endowed with 
a constant normal coefficient of restitution, e. Therefore 

y'n = ~^'"n, (20) 

where u„ and v'^ denote the normal relative velocities of two 
colliding particles before and after the collision respectively. 
The particles' diameter we denote by d, which sug- 
gests a dimensionless parameter analogous to g. For con- 
sistency with previous work (Ohtsuki 1993, Canup and Es- 
posito 1995), we use the ratio 






2Gm 
3dMV 



-1/3 



(21) 



where ru is the mutual Hill radius of the particles. The pa- 
rameter rp may be recast as 

, 1/3 

^2i/3 P£ 



(22) 



where pp is the mass density of a particle, pc is the mass 
density of the central body, with Re denoting its radius. In 
the case of Saturn's rings, it follows that r-p < 1 at radii 
greater than roughly 125,000 km, if we take pp to be that of 
crystalline water ice (900 g/m''). At the radius of the F-ring, 
rp falls to 0.89. 

Finally we can derive the ratio of the particle radius and 
initial particle spacing from d/h = 0.873 r-p g^'^ ■ Because we 
must have d < h (or else particles overlap), there is a re- 
striction on the size of rp given g, which is easy to calculate: 
rp < 1.1455"^/^ 



3.2 Clustering criteria 

The key to the long-term evolution of the ring is the ca- 
pacity of ring particles to form long-lived aggregates. In this 
subsection we discuss a number of simple criteria that might 
help us understand the simulation outcomes. 

Initially we consider whether simple aggregates of two 
particles in contact can resist tidal disruption. The simplest 
case is that of an aggregate in synchronous rotation, which 
appears as a non-rotating aggregate in the frame of reference 
used in this paper. The particle locations are given by t/i = 
2/2 = and xi — —d/2 and X2 ~ d/2. By testing the limiting 
case of no contact force, Eqs ([l]) and H tell us that this 
configuration is an equilibrium solution if 



s; 1, 



(23) 



(Weidenschilling et al. 1984, Longaretti 1989). 

This 'force' condition provides a clear criterion, but it 
applies only when two particles are touching in this special 
configuration. How particles ever reach such an arrangement 
is not addressed: they may not be all that likely in real 
coUisional systems. Consequently, force criteria may over- 
estimate clump formation. Another approach is taken by 
Ohtsuki (1993) and Canup and Esposito (1995) who exam- 
ine the energetics of a collision via the Jacobi integral and 
compute a condition for accretion during a binary collision. 
These criteria relate Vp and e. Canup and Esposito find that 
when 



e < 



\ + (2/3)r2 - 9 
-h c2 + (2/3)r2 ' 



(24) 



colliding particles should aggregate. Here Ve is escape ve- 
locity and c is the velocity dispersion of the system. Both 
are scaled by r^Q,, and so v^ = \/G/rp. If we assume that 
c ~ lie in equilibrium, then we have a relationship between 
e and r^ that tells us whether a typical collision results in 
aggregation, and hence whether our system has a tendency 
to become clumpy. In the limit of perfectly dissipative col- 
lisions (e = 0), which is the most conducive to aggregate 
formation, the criterion becomes 



rp s; 0.691, 



(25) 



signiflcantly less favourable than ( 23 1 



These criteria hold only for binary aggregates or bi- 
nary collisions; but if a binary (or larger) clump has already 
formed then additional particles can attach themselves to 
it with greater ease (Weidenschilling et al. 1984). It follows 



that if criterion ( 23 1 is satisfled and yet ( 24 1 is not satis- 



fled, a rare event may occur whereby a long-lived aggregate 
forms and then successfully absorbs other particles, thereby 
growing larger and larger. For fixed e, it follows that an inter- 
mediate regime may exist, with rp lying between the limits 
obtained from ([24| and ([23| (Salo 1995). For e = 0, this 
would be 0.691 < r^ < 1. In this regime long-lived clumps 
do not form in general — almost all collisions do not re- 
sult in capture — but given sufficient time a rare particle 
encounter occurs that does result in an aggregate. This ag- 
gregate may then grow successfully, an embedded blob in an 
otherwise monodisperse ring. 



3.3 Numerical method 

Using the coUisional N-body code REBOUND (Rein and 
Liu 2012) we study the gravitational instability by directly 
evolving Eqs. ([l])-(l3|. Because of the small number of par- 
ticles, we can use direct summation to accurately calculate 
self-gravity and do not have to use a tree or FFT-based 
gravity solver. 

We use the mixed-variable symplectic epicycle integra- 
tor SEI (Rein and Tremaine 2011) which is ideally suited 
for this study. A symplectic integrator does not introduce 
artificial trends in formally conserved quantities. The mixed- 
variable integrator gives a large accuracy gain when the par- 
ticle motion is dominated by epicyclic motion. Because of 
these properties, we have numerically converged results in 
most simulations using a time-step as large as one tenth of 
the dynamical time, dt — 10~^r2~^. 

Usually, one employs several ghost boxes in local N- 
body simulations to ensure that there is no special place 
in the shearing sheet. The more ghost boxes are used, the 
better. Here, we find it is sufficient to use only one ghost 
box in both the positive as well as the negative y direction. 
There are no ghost boxes in the x direction, since we are 
dealing with a string of particles. Thus each particle feels 
the attraction of 3N — 1 other particles. 

We pre-calculate the residual force error Ei in the calcu- 
lation of fy comparing the initial numerical setup to a truly 
infinite string of particles with perfect periodic spacing h, 
cf. Eq. H . This is similar to the idea of Ewald summation, 
and the error can be written down explicitly in terms of the 
first derivative of the digamma function 



E, 



Gm 



{ij)' {2N - i) ~ iP' {I + i + N)) 



(26) 



where i is the index of the particle running from (left) to 
N—1 (right), and a prime denotes differentiation. This error 
is then subtracted from the numerically calculated forces at 
every time-step. Using this trick, artificial effects are min- 
imised at the box boundaries even when the instability is fol- 
lowed over many dynamical timescales and grows by many 
orders of magnitude. 

We use an instantaneous collision model. Hence mul- 
tiple collisions during one time-step may not be treated 
correctly. Furthermore, we use a minimal impact velocity 
(0.05 dfl ) so as to avoid overlapping particles when aggre- 
gates form. This velocity scale is set much smaller than the 
velocity dispersion and so should not affect the outcome. 



3.4 Linear regime 

As a numerical check, the analytic growth rate given by 
Eq. ( [TT| ) is verified using the N-body code. We take A'' = 200 
particles initially lined up in a string and simulate their early 
evolution. The particles' location and velocity are perturbed 
with the unstable eigenfunction of fastest growth. The initial 
perturbations possess amplitudes of 10"^ and we stop the 
evolution when they have grown by 2 orders of magnitude. 
In Fig|4]we plot the largest growth rate as a function 
of the dimensionless parameter g. We also show the analytic 
result from linear theory. One can see that the agreement is 




Figure 4. Largest growth rate of the system in numerical simu- 
lations and linear theory as a function of the parameter g. Note 
that the detailed structure at small g is captured perfectly. Here 
N = 200. 



Case 


9 


rp 


Instability 


Aggregation 


d/h 


(i) 


0.1 


0.2 


Oscill. 


Likely 


0.081 


(ii) 


0.1 


0.8 


OsciU. 


Unlikely 


0.32 


(ill) 


1 


0.2 


Monot. 


Likely 


0.17 


(iv) 


1 


0.8 


Monot. 


Unlikely 


0.7 



Table 1. Parameter sets for the four runs discussed in detail in 
Section 3.5 



excellent, which provides a useful validation of the numerical 
integrator, and also a check on the analytic theory. 

3.5 Nonlinear evolution 

Long time simulations were ran for a variety of parameters, 
though first we concentrate on the following four parame- 
ter sets presented in Table I. In the table, 'OsciU. /Monot. ' 
refers to the form of the initial instability, oscillatory or 
monotonic (determined by g), and 'likely' or 'unlikely' ag- 
gregation refers to the likelihood of clumping as decided by 
(241). In all cases we set e = 0.5 and N = 200. Note that 



the 



0.8 cases fall within the 'intermediate regime' dis- 



cussed earlier: aggregates can still exist in principle. These 
parameter choices hence let us observe the potentially dif- 
ferent outcomes of the two forms of gravitational instability, 
while also testing the clustering criteria of Section 3.2. 

We find that on medium to long times the sinrulations 
yield two typical evolution tracks depending on the param- 
eter Vp. It appears that after a few orbits the ensuing dis- 
ordered state retains little memory of the linear instability 
that gave rise to it. We take simulations (i) and (ii) as ex- 
amples of these two tracks. 

First we treat case (i), as it is the simplest. These pa- 
rameters yield oscillatory instability and physical collisions 
that usually lead to capture. In Fig. [5] screenshots of the 
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Figure 5. Snapshots of a simulation of gravitational instability. Parameters are: g = 0.1, Vp = 0.2, e = 0.5, N = 200. Note that only a 
central portion of the computational domain is shown. Particle sizes have been inflated to aid their visualisation; hence particle aggregates 
appear as overlapping particles. The simulation is initialised by small amplitude random noise. 



evolution are presented corresponding to different times. In 
order to better visualise the results, these images only repre- 
sent the central portion of the computational domain, com- 
prising initially of 40 particles. Particles are represented by 
black circles and their radius has been inflated so they are 
more easily seen. In panel 1, the equilibrium set-up is shown, 
while in panel 2, the system is on the verge of departing from 
the linear regime of the oscillatory instability. However, the 
perturbations are still too small to be easily observed. Panel 
3 shows the outcome of the first round of collisions. These 
initial collisions induce accretion and precipitate disordered 
non-planar motion in the resulting aggregates. Gravitational 
encounters liberate orbital shear energy and the ring 'warms' 
up with a velocity dispersion ~ Ve', in addition, angular 
momentum is transported and the ring spreads, its width 
eventually oc t^ , for p close to 1/2. Particles and particle 
aggregates continue accreting through physical collisions, in 
accordance with the criterion (24 1, and after 8 orbits the 



number of bodies has decreased significantly. Because we 



have infiated the particle radiuses in the images, the aggre- 
gates appear as overlapping circles. 

The evolution of case (ii) is represented by eight screen- 
shots in Fig.|6] Unlike case (i), the system does not steadily 
collapse into a smaller set of aggregates. Instead it swiftly 
degenerates into a disordered spreading state characterised 
by a velocity dispersion, as before, of order v^ (consistent 
with Ohtsuki 1999 and Ohtsuki and Emori 2000). The eight 
panels clearly describe this spreading. However, what they 
also show is the continuous generation and dissolution of 
small particle aggregates. Particles come together and then 
are disrupted by tidal forces, spin, or forceful collisions with 
other particles. These temporary aggregates resemble the 
dynamic ephemeral bodies (DEBs) proposed by Weiden- 
schilling et al. (1984). Alternatively, we may think of them 
as analogues of the gravity wakes exhibited by optically 
thicker rings (Salo 1992). We have found that some aggre- 
gates persist for the length of the simulation. These bodies 
may have become sufficiently large to pass the tidal disrup- 
tion threshold for these parameters. Simulations were also 
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Figure 6. As in Fig. [5] but with parameters g = 0.1, r^ = 0.8, e = 0.5, N = 200. 



run for rp = 0.85, 0.9, and 1 and long-lived clusters failed 
to appear, though occasional two-particle aggregates were 
spotted. 

These runs confirm an intermediate regime of rp be- 
tween continuous aggregation and no aggregation at all. This 
regime is characterised, in particular, by temporary aggre- 
gates and by the potential formation of persistent growing 
clusters that resist disruption amidst the surrounding melee. 



3.6 Clumping criterion 

In order to measure the efficiency of clustering in detail we 
conducted a parameter survey of e and Vp. We employ a 
convenient 'one-dimensional' measure of the system's degree 
of clustering which we now describe. Particles are sorted 
along the azimuthal direction and the azimuthal distances 
between neighbouring particles Ay calculated. (Their rela- 
tive radial distances are discarded.) The mean of the result- 
ing distance distribution (Ay) stays constant throughout a 
simulation and is equal to h. However, its standard devia- 
tion will evolve over time. We average this standard devia- 
tion over the length of the run and scale it by h; the final 



averaged quantity we denote by q and associate it with the 
degree of aggregation throughout the run. Note that g = 
in the perfectly ordered equilibrium state that we start with, 
while g = 1 if the particles are distributed entirely randomly. 
In the latter case, Ay exhibits an exponential distribution 
with scale parameter h. Aggregation corresponds to g > 1, 
with q taking its maximum value q ~ \/iV when all the 
particles lie in a single cluster. In our simulations, we find 
that q never exceeds about 6. The one-dimensional measure 
q is well-suited to a narrow spreading ring, and is also sim- 
pler to implement than the two-dimensional Renyi entropy 
employed in Karjalainen and Salo (2004). Like the Renyi 
entropy q cannot distinguish between the continuous forma- 
tion of many temporary clusters and permanent aggregates. 
Figure [7] plots q in greyscale as a function of e and Vp. 
It summarises 2240 separate simulations each run up to 200 
orbits. The figure illustrates clearly the hard boundary at 
rp = 1; for larger rp, we find q — 1 and the particles are 
distributed randomly with little aggregation. Moreover, our 
results show that aggregation is insignificant until at least 
Vp ~ 0.8, though this limit decreases for larger e. The q value 
jumps abruptly near this critical rp, from slightly greater 
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Figure 7. The measure of dumpiness g as a function of e and rp. 
Superimposed on the greyscale plot is the Canup and Esposito 
criterion (1995) for gravitational capture after a binary colUsion 
(the black line). The other parameters are g = 0.1 and TV = 200. 



than 1 (little accumulation), to g ~ 5, which we associate 
with the onset of continuous accretion. In addition we have 
overplotted the Canup and Esposito criterion (1995) given 
in Eq. ( |24[ ), which appears to underestimate the prevalence 
of clustering, especially for higher e. 

At small Tp it appears that q decreases slightly. This 
is perhaps due to the fact that in this regime the particles 
are so small that they rarely collide with each other. Ex- 
tremely long time integrations are hence required to reach 
the systems' accretion time-scales. Our simulations, which 
were run up to 200 orbits, were perhaps too short. It is also 
possible that accretion is suppressed in hot rings composed 
of very small particles. 

The numerical results are compatible with self- 
gravitating simulations of broad rings (Salo 1995, Kar- 
jalainen and Salo 2004), which also generate long-lived ag- 
gregates. However, the optical depth is much greater in those 
simulations, and they are significantly perturbed by gravity 
wakes — conditions that make aggregate formation more 
likely and more rapid than in our work. In particular, Kar- 
jalainen and Salo (2004) find that when r = 0.25 continuous 
accretion occurs for rp < 0.87 if e = 0.1, and for rp < 0.85 
if e = 0.5. These critical values are larger than ours, but at 
lower r their critical rp approaches 0.8, which is closer to 
our critical value. It should be noted that these low r runs 
were undertaken with a varying coefficient of restitution. 



4 CONCLUSION 

In this paper we have studied the gravitational and coUi- 
sional dynamics of a stream of co-orbital particles, plotting 
the course of the system's evolution from initial instability 
to its nonlinear saturation. Gravitational instability attacks 
the equilibrium in one of two ways: as either a monotonic 
clumping or as growing epicyclic motions. The latter's mech- 
anism is particularly interesting, as it consists of the inter- 



action of gravitational attraction and angular momentum 
exchange. In this way it shares some properties with other 
disc instabilities, such as the Papaloizou-Pringle instability 
and the MRI. 

The nonlinear coUisional evolution of the instabilities 
are tracked with N-body simulations. These show that the 
final saturated state is relatively insensitive to the form of 
the initial instability and is instead controlled by the inelas- 
ticity of collisions (measured by e, the coefficient of resti- 
tution) and the ratio of particle diameter to the mutual 
Hill radius (rp). The system tends to move towards three 
regimes: (a) for r^ > 0.83 a hot disordered fiow ensues with 
little or no clumping; (b) for rp less than 0.83 but greater 
than a second critical value rc{e), the system continuously 
forms short-lived small aggregates and occasionally a per- 
manent cluster; (c) when rp < re most collision result in 
accretion and the system accumulates into permanent and 
continuously growing conglomerates. 

Regime (b) is especially relevant to conditions in the 
F-ring of Saturn, where the rp parameter is close to 0.83. 
Our simulations, like Salo (1995) and Karjalainen and Salo 
(2004), are hence consistent with observations of embedded 
large bodies in the F-ring, and theories that attribute their 
prevalence to a size distribution dynamics comprising grav- 
itational aggregation, and tidal and coUisional disruption 
(Barbara and Esposito 2002, Esposito et al. 2011). Our sim- 
ulations directly animate these processes from first princi- 
ples, though the real system exhibits many more physical 
effects than our model. These include 'stirring' by moons, 
adhesive surface forces, and a size distribution. 

Finally, the linear instabilities presented here may have 
analogues in dense narrow rings that exhibit orbital shear 
(Goodman and Narayan 1988, Papaloizou and Lin 1989). 
Gravitationally unstable streams of material torn from 
tidally disrupted satellites may be the progenitors of these 
narrow rings, including the Uranian e-ring (Leinhardt et 
al. 2012). Allowing for the modifications wrought by shear, 
the basic mechanism of instability should be much the same, 
combining gravitational clumping and angular momentum 
exchange. The simple system we study here, which permits 
us to isolate and understand the salient physics, can then 
provide an inroad into the more complicated dynamics of 
the confined shearing system. 
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